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ABSTRACT 

We present a linear analysis of inviscid, incompressible, magnetohydrodynamic (MHD) shallow 
water systems. In spherical geometry, a generic property of such systems is the existence of five wave 
modes. Three of them (two magneto-Poincare modes and one magneto-Rossby mode) are previously 
known. The other two wave modes are strongly influenced by the magnetic field and rotation, and 
have substantially lower angular frequencies; as such, we term them "magnetostrophic modes" . We 
obtain analytical functions for the velocity, height and magnetic field perturbations in the limit that 
the magnitude of the MHD analogue of Lamb's parameter is large. On a sphere, the magnetostrophic 
modes reside near the poles, while the other modes are equatorially confined. Magnetostrophic modes 
may be an ingredient in explaining the frequency drifts observed in Type I X-ray bursts from neutron 
stars. 

Subject headings: waves — magnetohydrodynamics — stars: neutron — X-rays: bursts 



X 
o 



> 

On 
p 

00 
O 
On 
O 



1. INTRODUCTION 

Shallow water wave systems are those in which the 
transverse length scales considered are much larger than 
the water height. They have been extensively studied, 
have a wide range of physical applications, and have 
been used to understand terrestrial and planetary sys- 
tems (Matsuno [1966]; Gill [1982], hereafter G82; Ped- 
losky [1987]; Braginsky [1998]; Holton [2004]; and Kundu 
& Cohen [2004], hereafter KC04). Shallow water equa- 
tions capture the large scale dynamics of thinly-stratified 
atmospheres under the influence of rotation and thermal 
forcing. 

In astrophysical settings, shallow water models were 
used to study the spread of accreted matter onto neu- 
tron stars (Inogamov & Sunyaev 1999), and the ignition 
and propagation of Type I X-ray bursts in neutron star 
atmospheres (Spitkovsky, Levin & Ushomirsky 2002). In 
addition to fast rotation and strong gravity, many neu- 
tron stars in accreting systems also possess appreciable 
(< lO^G) magnetic fields. Tension in the fields induced 
by horizontal motions in the atmosphere can affect the 
dynamics, and thus it is interesting to study the effect 
of frozen-in magnetic fields on the behavior of the atmo- 
sphere. 

Gilman (2000) pioneered the use of magnetohydrody- 
namic (MHD) shallow water systems for studying the 
solar tachocline, which inspired several follow-up studies 
(e.g., Schecter, Boyd & Gilman 2001; Zaqarashvili et al. 
2007, 2009). Most works on MHD shallow water systems 
consider initially toroidal magnetic fields and/or slow ro- 
tators, as is applicable for tachocline research. In neu- 
tron star-related applications, however, it makes sense 
to consider initially vertical (or radial) magnetic fields in 
the presence of fast rotation. We study this case in the 
present paper. 

The fundamental governing equation for shallow water 
waves on a sphere is known as Laplace's tidal equation. 
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Fig. 1. — Qualitative behaviour of the wave modes, as controlled 
by the MHD Lamb's parameter. The governing equation for the 
latitudinal velocity perturbation, in the respective limits, is given 
in parentheses. 



By linearizing this equation, one can obtain solutions 
to the water height and velocity perturbations, as well 
as dispersion relations for the angular frequencies of the 
waves. In rotating systems, it was realized (Longuet- 
Higgins 1965) that an important quantity in such studies 
is Lamb's parameter. 
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(1) 



where TZ is the Rossby number'^. Physically, shallow wa- 
ter systems with large values of e are "fast rotators" . In 
a seminal paper, Longuet-Higgins (1968, hereafter LH68) 
explored solutions to Laplace's tidal equation for a wide 
range of values for e and demonstrated that analytical 
forms exist for e ^ and |e| — *■ oo. While e < so- 
lutions may appear unphysical, LH68 realized that they 
are relevant to the study of forced oscillations. 
When magnetic fields are considered, one needs to in- 

Defined as the ratio of inertia to Coriolis forces on the scale of 
the planet or star. 
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stead examine the MHD Lamb's parameter, 

i^I{B,,n)e. (2) 

The fmiction T — T{Bz,^), which we wiU define later, 
depends on the rotation frequency of the system Q and 
the magnetic field strength B^; for hydrodynamic waves, 
T = 1. Unlike in the case of hydrodynamic systems, e < 
waves are directly relevant in MHD — for — e ^ 1, they 
reside near the poles on a spherical surface. By contrast, 
e ^ 1 waves reside near the equator. Figure [T] shows a 
schematic demarcating the asymptotic behaviour of the 
wave solutions and the limiting equation that governs the 
latitudinal velocity perturbation. 

In this paper, our goal is to demonstrate that the 
methods of LH68 can be generalized to obtain analyt- 
ical, asymptotic (|e| ^ 1) solutions to the MHD ana- 
logue of Laplace's tidal equation. These solutions can 
then be used as a guide towards obtaining |e| ~ 1 solu- 
tions, which must be computed numerically. Throughout 
the study, we shall adopt a radial magnetic field for sim- 
plicity (and algebraic amenability). Our results can be 
straightforwardly generalized for arbitrary magnetic field 
configurations. 

Readers unfamiliar with the classical shallow water 
analysis are referred to Appendices [X] and |B1 where the 
hydrodynamic treatment is presented in Cartesian and 
spherical geometry, respectively. In fj^l we add magnetic 
fields to the Cartesian analysis. Our efforts culminate in 
fJJl where we explore MHD solutions on a sphere. We 
discuss the implications of our results in Sj4l A concise 
summary of the paper is presented in SJS] 

2. MHD SHALLOW WATER WAVES: CARTESIAN 
COORDINATES 

Consider a Cartesian coordinate system in which a 
shallow layer of water of height h — h{x,y,t) resides 
at z > on an infinite plane. The fluid moves with a 
velocity v — [v^^ Vy,Vz) such that 



magnetic field: 



Vx = Vx{x,y,t), 
Vy = Vy{x,y,t). 

The fluid is assumed to be inviscid and incompressible 



(3) 



2.1. f -Plane Treatment (Magneto-Poincare and 
Geostrophic Waves) 

2.1.1. Equations 

Consider the Euler equation with the Coriolis force and 
magnetic tension terms (e.g., Draine 1986): 

1 



Ov 1 

— + v.Wv = VP -2{nxv 

at p 



Incompressibility implies (see Appendix lA.ip 

|^ + V.(..) = 0. 

The induction equation, in the ideal MHD limit, provides 
an additional equation of motion: 



(5) 



OB ^ ^ 

— = V X [v-x. B 
dt \ 



(6) 



We consider small perturbations to the water height, 
as well as to the x- and y-components of the velocity and 



Vx 

K 



h = H 
Bx = Bxo 
Bi, — B. 



h' 



(7) 







(8) 



The field is initially vertical (we will set Bx^ = By 
later), and is frozen at the bottom of the atmosphere. 
The largest gradient of the field is due to the vertical 
shear in the layer, which we approximate as: 

dB^^ B^ 
dz H ' 

where the index j represents the set j = {x,y}; B^ is 
taken to be constant. The net effect of these approxi- 
mations is to introduce restoring tension forces that pull 
the fluid elements back into their original horizontal po- 
sition. In the shallow water formalism, one deals with 
height-averaged horizontal velocities and accelerations. 
Collectively, it is then reasonable to take the restoring 
force (per unit volume) to have a mag nitude B'^/AttH 
and a direction opposite to that of the horizontal dis- 
placement. 
The following wave solutions are sought: 



exp i {kxX -f kyy 



y 
h' 



Vyg expi{kxX 
ho exp i {kxX 4 
bx = bxo exp i [kxX - 
by = byo exp i {kxX - 



~ujt), 

-Lot), 
yy - Lot), 

kyy-ujt), 
kyy-ujt), 



- kyy 



(9) 



where the wave vector is fc = (kx,ky,0). Linearization 
of the equations of motion yields a 5 x 5 matrix (see 
Appendix lA.ip : 
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A 



where 
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An 



Ai = 
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Bz 



kxVx 
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2± 



Byo^y) , 



(10) 



Ai±- {Bxakx 

Airp 

A3 = AnpA2+ + iBxgkx, 
Ai = ATTpA2+ + iBygky. 

A practical note about evaluating detA is that one is free 
to permute the rows of A, and they should be arranged in 
a way so as to minimize the number of sub-determinant 
(i.e., of 4 X 4 matrices) evaluations. 

rsion Relation 



2.1.2. 

Evaluating detA — 0, setting Vx — Vy = Bx^ = Byg 
and collecting terms yields the dispersion relation: 



Lo* - uj^ [gHk^ + + 2{vaIH)'^] 
+ {VA/Hf [gHk^ + (va/H)^] = 0, 



(11) 
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where va = i?^/2y^7rp is the Alfven speed. Its solution 
is 



(12) 

Some intuition can be developed by examining the so- 
lutions in 17 = limit: 



2 _ {gHe + {vAlH)\ 



\{VA/Hf 



(13) 



The first mode is the longitudinal "magnetogravity 
wave" (Schecter, Boyd & Oilman 2001), where the restor- 
ing force is a combination of the rising fluid height and 
magnetic pressure increase due to the compressive mo- 
tion of the gravity wave. The second mode, which is 
non-dispersive, represents the Alfven wave with no height 
perturbation and torsional oscillation in the transverse 
direction, whose restoring force is magnetic tension. 

When we revert to rotating systems, both modes 
begin to couple longitudinal and transverse velocities 
through the Coriolis force, but have different polar- 
izations. The higher-frequency mode in equation \12\ 
is the rotationally-modified magnetogravity wave, or 
"magneto-Poincare" wave — the magnetic tension term 
in equation ([4]) adds in phase with the Coriolis force, 
thus enhancing the restoring force and speeding up the 
oscillation. In the lower-frequency mode from equation 
p2|) . magnetic tension tries to balance the Coriolis force. 
To emphasize this balance, we call this the "magne- 
tostrophic mode" . 

Both magneto-Poincare and magnetostrophic modes 
come in eastward- and westward-propagating vari- 
eties, giving the four modes as required by equation 
pT|) . The eastward- and westward-propagating waves 
of each branch have the same frequency by absolute 
value. As the magnetic field is turned off, the magne- 
tostrophic mode disappears, while the Poincare mode (or 
rotationally-modified gravity wave) survives as expected. 



2.1.3. Length Scales 

In hydrodynamic systems, rotational effects become 
important at wavelengths of A = 1/fc > Ar, where 

|fc| = 



k 



1,2 



and 



Ar = 



2n 



(14) 



is the "Rossby radius of deformation" (e.g., chapter 7 of 
G82). The Rossby radius can also be obtained by arguing 
that there exists a radius at which the radial fluid flow 
is diverted by the Coriolis force, i.e., X^^l ^ co, where 



Co 



9H 



(15) 



is the shallow water wave speed. Physically, an adjust- 
ment to "geostrophic balance" occurs at the Rossby ra- 
dius on a timescale ^ l/fl (pg. 201 of G82). 

In MHD systems, even without explicitly calculating 
the group velocity, one can see from equation (jlip that 



if Va ~ flH, we have 



^RB 



9H 



4f]2 + 2{va/H) 



(16) 



The quantity Arb can be understood in the following 
manner: if the forces due to rotation and magnetic 
tension are equally important, then at wavelengths of 
A ~ Arb they balance out the effect of gravity. 

A somewhat more relevant quantity to define is the 
Alfven radius, Ab. In the absence of rotation, the forces 
due to magnetic tension and gravity balance at A ^ Ab, 
where 



Ab = Arb {n = 0) 



H. 



(17) 



The preceding expression can be approximately obtained 
by arguing that Ab ~ cqH/va- The relative importance 
of rotation and magnetic tension can be judged from the 
ratio Ar/Ab. 

2.2. P -Plane Treatment (Magneto-Poincare and 
Magneto-Rossby Waves) 

By analogy with hydrodynamic systems, we term the 
very slow waves with large wavelength "magneto- Rossby 
waves" . Generalizing the analysis in Appendix IA.31 we 
set Vx = Vy = Bxo — Byg = and consider the following 
set of linearized equations: 

d< 
dt 



dt 
dh' 
'dt 
dbx 

dt 

dby 

dt 



= 0, 



9< 
dx 

H 
H 



dy 



Alby = 0, 



(18) 



where f — 2^1 sin 8 -t- /??/ is the Coriolis parameter and Q 
denotes the latitude. We differentiate the first equation 
in (|18p with respect to y and seek wave solutions from 
the entire set of equations^. We keep only first order 
terms in the Coriolis parameter, eliminate bxo and byg, 

and construct the A matrix: 



A = 



I 



V 



K H , 



(/3 + %kyh) 



h 

KB 



— — LO 



kyH 



gkx ky 
igky 



\ 



J 



(19) 



As before, setting detA — yields the dispersion relation, 

uj'^ - uj^ \gBk^ + /o + 2 {va/B) 

- gBkxf3u; + [va/B^ [gBk^ + (va/B)''] = 0, 

which reduces to equation (fTTj) when /3 — as expected, 
i.e., the /9-plane approximation reduces to the /-plane 

* As a check, we are able to reproduce equation (15) of Za- 
qarashvili ct al. (2007) using our approach. 
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one at the poles. For very slow waves (uj <C /) with large 
wavelengths {kH <C 1), i.e., magneto- Rossby waves, the 
phase speed is 



where R is the radius of the sphere and 6 — 90° — is 
the co-latitude. We define the following quantities^: 

fi = cos 9, 



CR : 



where 




Ci = /3A^/C2, 

C2 = sin2e + (AR/AB)' 

Xj: = 1/kx- 



(20) 



(21) 



D = -sin( 

_ LU 

^ = — 



dd 



V4, = v'^ sin 9, 
bg = be sin 6*, 
L = brh sin6. 



(24) 



The negative sign implies that the phase propagation is 
westward. The root of the dispersion relation with the 
slowest angular frequency has the phase speed: 




Ci 



(22) 



The introduction of vg and allows one to avoid singu- 
larities associated with 9 = 0°. 
We seek the wave solutions. 



ve — vgg exp i{s4> — ix)t), 
V4> = exp i {s(t> - ujt), 
h' = exp i (s4> — Lot) , 



(25) 



and eliminate the amplitudes for bg and b^ to obtain: 



When B2 ^ (Ab 00), we have Cp 



0. When 



P = and Ar /Ab ^ 1, we have Cp^ <C cq. This root can 
plausibly be associated with the east magnetostrophic 
mode. However, when considering the equatorially- 
confined, non-planar eigenfunctions on the /3-plane (see 
Appendix [Djl. we do not find this mode, yet do recover 
the two magneto-Poincare and magneto-Rossby modes. 
This suggests that the magnetostrophic modes do not 
reside near the equator on the sphere, and the /3-plane 
treatment may be inconclusive for these modes. We will 
confirm this in the next section. 



3. MHD SHALLOW WATER WAVES: SPHERICAL 
COORDINATES 

3.1. Equations 

By analogy with equation ([5]) , we specify a radial mag- 
netic field and allow for perturbations in the 9— and (f>- 
directions (bg and b^, respectively). The linearized equa- 
tions of motion become: 



dig a : 9 ■ a^h' Brbg 

— -2f7.,cos0-F-sm^— 

dv^ g dh' BrKp 

— -H217.,cose + -— + ^^0, 

' '~dt^ R '89 ^~R^ 

dbg _ VgBr 

'dt ~ ~ir' 

db^ _ v^Br 
~dt ~ H ' 



0, 



0, 



(23) 



Xvg„ + fiv^g + D1J0 = 0, 

fJ-VOo + XV - Sl]o = 0, 

zue (1 - ^^) 770 - Dvg„ - sv^„ = 0, 



(26) 



where s is the toroidal wavenumber, vg^ = ivg^ and 
770 = gho/2nR. Notice that the spherical hydrodynamic 
amplitude equations [equation (|B2p ] and equation (j26|) 
have identical structures, except that 



(27) 



takes the place of n7 in a couple of places, with va 
Br/2y/Trp being the Alfven speed and 



■UJA 



VA 



Ar 

2nH V2Ae 



(28) 



The equation for vg„, obtained from the set (|26p . has 
the form 



LsVg„ = 0, 



(29) 



where the operator is 



SjJL 



(30) 



In the place of e, we now have e = vot/x in the MHD 
case. 

® One of the earliest papers to use at least some of these trans- 
formations is Margules (1893). 
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3.2. Near-Equator Solutions (e 3> 1 j 

When e 3> 1, the governing equation for vg^ reduces 
to the spheroidal wave equation with the separation con- 
stant, 



As„ (q) = eX , 

X 



where q = vl. The dispersion relation is 



(31) 



w 



\ e 



vo\ {21 + 1) 



137 



ZD -\- W A = 0. 



(32) 



Eight roots exist for equation ([5^ , but only three of these 
correspond to the angular frequencies for the magneto- 
Rossby and magneto-Poincare modes. We checked this 
by requiring the eigenfrequencies and their corresponding 
eigenfunctions, which will be derived shortly, to satisfy 
the set of equations in in the appropriate limits. 
We define 

( (33) 
^ = exp I — — I exp i {s4> — u)t). 



Notice that as the magnetic field strength increases, the 
waves are "pinched" closer to the equator because of the 
the exp {-X'^/2) term. 

Near the equator (|/j,| <C 1), the spheroidal wave equa- 
tion reduces to the parabolic cylinder equation. Its solu- 
tion is 



-ifia Hi 



(34) 



where Hi = Tii{X), 'Hi{X) is the Hermite polynomial, 
I = n — s and n plays the role of the poloidal wave num- 
ber. The quantity 770 = gho/2D,R is determined once 
the normalization Hq is specified. Using the amplitude 
equations in (j26p . we have: 



^0 



ex^ 



(35) 



Also, 



iHi-i — -^i+i 



(36) 



modes share the same eigenfunctions: 

, _ sf]o 



h' 



Br 



srjo 



ujH (ex^ - s^)x 
ho 



eX^ — 



l'Hi-iA+ + -Hi+iA- 



where we have 



(37) 



(38) 



The east magneto-Poincare mode is described by the 
solutions: 



bif, 
h' 



ho 



e-X^ - s 



2 y IHi-iA- + -Hi+iA^ 



(39) 

It is worth noting that the equatorial eigenfunctions 
can be obtained in the /3-plane approximation as well 
(see Appendix iPjl . 

3.3. Near-Pole Solutions (^e 3> 1 j 

Near the poles, the governing equation for the latitu- 
dinal velocity perturbation is Whittaker's equation (see 

Appendix IB.Sp . The Lg operator in the hydrodynamic 
[equation (|B3|) ] and MHD [equation ([50)1 ] cases are math- 
ematically identical, except that e and x take the places 
of e and tu, respectively. Therefore, the analysis per- 
formed in Appendix IB . 31 can be identically applied here. 
However, the ^ ~ — 1 solution yields the east magne- 
tostrophic mode, because 



X 



-1 



TO « - -1 
2 



1 + 4ti7^ 



(40) 



Likewise, the x ~ +1 solution yields the west magne- 
tostrophic mode. We assume the form 



X = ±l 



(41) 



from which the dispersion relation for the slow modes is 

— 2to^ j 

(42) 

it 2w j^w -\- I ujj^ 



0, 



The west magneto-Poincare and the magneto- Rossby where Q = 2m -f 2Z + 1 ; for x ~ ± 1 , we have m ~ | s ± 1 1 
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Fig. 2. — Solutions to the dispersion relations for MHD shallow 
water wave systems in spherical geometry (s = 1 and n = 2). 
Note that for the magneto-Poincarc modes, the curves for the three 
different values of Aj^/Ab = 0.01,0.05 and 0.1 overlap. 

The solution to equation (|29|) is 

ve„ « Vo exp (-^) f (f ) , (43) 



where Vq is an arbitrary normahzation constant, C 
is the associated Laguerre polynomial and 

The other amplitudes can be computed using: 



(m; 



6^ (1 



(44) 
(45) 



X 



Knowledge of the amplitudes allows one to compute the 
magnetic field perturbations by taking the real parts of 
the following expressions: 



b0 = 



iBrV'a 



(46) 



Only two of the roots in equation (|42)) are physical; we 
checked this by again requiring the eigenfunctions to sat- 
isfy the set of equations in (j23p in the appropriate limits. 

4. DISCUSSION 

4.1. Existence of Wave Modes 

Our analyses in the previous sections have shown that 
five wave modes exist in MHD shallow water systems. 
We now discuss a more intuitive way of understanding 
why they exist. Firstly, consider fluid flow on the surface 
of a non-rotating, non-magnetized cylinder with gravity. 
Only gravity waves exist and the flow can generally be 
eastward- or westward-propagating. The Rossby mode 
does not exist even if we rotate the cylinder (hence pro- 
ducing the Poincare modes), because it requires the pres- 
ence of a latitudinally-varying Coriolis force. 



If we now allow for the presence of a magnetic field, the 
magnetostrophic modes appear alongside the magneto- 
Poincare modes. A magnetized, self-gravitating cylinder 
— regardless of whether it is rotating — possesses four 
shallow water wave modes. If we replace the cylinder 
with a sphere, the magneto-Rossby mode appears. Shal- 
low water systems on a non-magnetized sphere only have 
three modes: the Poincare and Rossby modes. 

Note that unlike on the /3-plane, there are both 
eastward- and westward-propagating magnetostrophic 
waves on a magnetized sphere. As these waves are more 
concentrated towards the poles, they are less suscepti- 
ble to the effects of varying the Coriolis parameter over 
the sphere. Furthermore — unlike in the case of the hy- 
drodynamic Rossby mode — the magneto-Rossby mode 
exists even for /3 — > 0, where it then reduces to the west 
magnetostrophic mode. 

4.2. Angular Frequencies 

The dispersion relations for the MHD modes on 
a sphere, described by equations and (|^ . are 

parametrized by the ratio of Rossby (Ar) to Alfven (Ab) 
radii, 

Ar 
Ab 



0.02B8 



-l/2rr-lQ-l 



(47) 



where for illustration, we have adopted Bs = Br /lO^ G, 
pe = p/lO^ g cra-^, H3 = H/IO^ cm, ila = rad 
s~^, g = 2 X 10^^ cm s~^ and R = 10 km, following 
Spitkovsky, Levin & Ushomirsky (2002) in the context 
of neutron stars. For these values, we have Ar « 2 km 
and Ab ~ 100 km (w^ « 0.3 km s~^). 

In Fig. [21 we calculate the angular frequencies of the 
wave modes for various values of Ar/Ab- Three fea- 
tures are apparent. Firstly, the magneto-Rossby mode 
has an angular frequency intermediate between those of 
the magnetostrophic and magneto-Poincare modes. The 
magneto-Poincare modes are insensitive to Ar/Ab- How- 
ever, the magnetostrophic modes are sensitive to Ar/Ab, 
because {tul ~ zu\. Since zua oc Br, such a property 
can conceivably be used as a magnetometer, provided 
the other quantities are known. 

4.3. Visualization of the Wave Modes 

In Figs. OandlH we show two-dimensional plots of the 
MHD wave solutions for the magneto-Poincare, magneto- 
Rossby and magnetostrophic modes. We adopt s — 1 and 
n = 2 for illustration, but note that the toroidal (s) and 
poloidal (n) wave numbers have geometrical interpreta- 
tions, corresponding to n — |s| -I- 1 and 2|s| wave nodes 
(i.e., "zero crossings") in the latitudinal and longitudinal 
directions, respectively. 

The qualitative behaviour of the modes agrees with 
the description of Matsuno (1966). In particular, the 
magneto-Rossby mode can be understood in the follow- 
ing manner: the height gradient — and therefore the 
pressure gradient — is positive from the low pressure re- 
gion in the middle {(j> = 180°) to the high pressure region 
on the right (Fig. [5]). Since v'y cx dh'/dx, we get v'y > 
and fluid flowing away from the equator. The reverse 
argument holds for the height gradient (because dh' > 
and dx < 0) between the low pressure region and the 
high pressure region on the left {</) — 0°), resulting in 
v'y < and fluid flowing towards the equator. 
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west magneta — Polncare 






Fig. 3. — Wave solutions (s = 1 and n = 2) for the magneto- 
Poincare and magneto-Rossby modes. For illustration, we have 
adopted Br = 10* G, H = 10^ cm, p = 10^ g cm'^, O = 10^ 
rad s~^, g = 2 X 10^'* cm s^'^ and R = 10 km. The contours 
represent the ocean height with light and dark colours correspond- 
ing to positive and negative perturbations, respectively; the height 
normalization is arbitrary. The arrows represent the velocity field. 

A feature that is absent from hydrodynamic systems 
is the magnetic pinching of the waves. Increasing 
concentrates the magneto-Rossby and magneto-Poincare 
waves closer to the equator; the magnetostrophic waves 
are concentrated closer to the poles. Although our solu- 
tions formally hold only in the large magnetic field limit, 
we verified the concentration of eigenfunctions and the 
values of eigenfrequencies for different modes by direct 
numerical integration of the spherical equations in (|23p . 
We used the method of Ivanov (2007) to find the regu- 




FlG. 4. — Same as Fig. \3\ but for the magnetostrophic modes 
with Br = 10^ G. 

lar solutions at the poles by series expansion and then 
applied the shooting method to determine the eigenfre- 
quencies of the odd and even modes. 

4.4. Type I X-ray Bursts from Neutron Stars 

Type I X-ray bursts are non-catastrophic, thermonu- 
clear explosions occurring on accreting neutron stars in 
low-mass X-ray binaries. They have typical rise and de- 
cay times ^ 1 s and ^ 10 — 100 s, respectively (e.g., 
Strohmayer & Bildsten 2006), and are often accompa- 
nied by millisecond (~ 300 — 600 Hz) oscillations (see van 
der Klis [2000] for a review) . As the X-ray burst reaches 
its peak luminosity and fades off, oscillations with fre- 
quency drifts ^ 1 to 5 Hz, relative to some asymptotic 
frequency, persist for a time 10 s. The asymptotic fre- 
quency is usually identified with the rotational frequency 
of the neutron star (e.g., Strohmayer et al. 1998). 

A popular interpretation is that these oscillations are 
caused by the rotational modulation of a growing "hot 
spot" on the surface of the neutron star. The timescale 
for thermonuclear burning is much shorter than the time 
needed to accumulate enough fuel to trigger thermal in- 
stabilities, allowing localized ignitions on the stellar sur- 
face. Initially, the hot spot moves backwards against 
the sense of rotation due to the (approximate) conser- 
vation of potential vorticity, much like westward- moving 
storms on Earth; it then spreads and engulfs the entire 
star (Spitkovsky, Levin & Ushomirsky 2002). Transverse 
pressure gradients cannot be ignored because the ignited 
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ocean of material is not fully degenerate (Gumming & 
Bildsten 2000), implying that an important characteris- 
tic velocity is the shallow water wave speed. 

Previous studies (e.g., Heyl 2004; Berkhout & Levin 
2008) have dealt with purely hydrodynamic shallow wa- 
ter systems and have focused on the Rossby mode ("r- 
mode"). It has been proposed that the frequency drift 
is caused by the Rossby mode, but it was quickly real- 
ized that the computed drifts are somewhat larger than 
what is observed. Our results in Fig. [5] are consistent 
with this view: for Ar/Ab ^ 0.01, the drift frequency is 
about V — —uj/2tt ^ 100 Hz (adopting the Spitkovsky, 
Levin & Ushomirsky [2002] parameter values). By con- 
trast, the west magnetostrophic mode has the frequency 
(in the limit that e ^ 1) 



47r 



Ar 

Ab 



Bl 



(48) 



which is of the correct order of magnitude (i.e., ~ 1 
Hz) if Ar/Ab 0.1 and Q ~ 10'^ rad s~^. However, 
if we assume Br and p to be constant with time, then as 
the accreted material cools and H decreases, v increases, 
contrary to what is required for explaining the temporal 
behaviour of the frequency drifts. Clearly, a full explana- 
tion of the frequency drifts also involves understanding 
the physics of accreting material onto the surface of a 
neutron star (Inogamov & Sunyaev 1999). 

It is not the intention of our present study to com- 
pute detailed models for the frequency drift. Rather, it 
is to point out that the west magnetostrophic mode is 
a promising candidate that deserves further attention. 
Note that the analytical eigenfunctions derived in 
are good approximations only for Br > 10^ G (with 
H ~ 10^ — 10"^ cm), whereas Type I X-ray bursters are 
believed to be ~ 10^ — 10* G systems. 

4.5. Future Work 

Our study of MHD shallow water waves is generic 
enough that it is readily amenable to various improve- 
ments. Natural extensions will be to consider viscous 
fluids with more general magnetic field configurations. 



The linear analysis can be generalized to a non-linear one 
via study of the Korteweg-deVries equation (e.g., KG04). 

5. SUMMARY 

The salient points of our study are summarized as fol- 
lows: 

1. We have performed a linear analysis of inviscid, 
incompressible MHD waves in the shallow water 
approximation. A generic feature of such systems 
is the existence of five wave modes: east and west 
magneto-Poincare; east and west magnetostrophic; 
and magneto-Rossby. 

2. Analytical functions for the velocity, height and 
magnetic field perturbations are obtained in the 
limit that |e| ^ 1. These functions arc useful as 
a guide towards obtaining numerical solutions or 
performing simulations. 

3. For reasonable values of the parameters, the 
magneto-Rossby and magneto-Poincare modes be- 
long to the e ^ 1 family of solutions and reside near 
the equator. The magnetostrophic modes belong to 
the — e 3> 1 family of solutions and reside near the 
poles. Increasing the magnetic field strength fur- 
ther concentrates the waves near the equator and 
poles, respectively. 

4. The west magnetostrophic mode is a potential can- 
didate for explaining the observed frequency drifts 
in Type I X-ray bursts from neutron stars. 
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APPENDIX 

A. CLASSICAL SHALLOW WATER WAVES: CARTESIAN COORDINATES 

A.l. Without Rotation (Gravity Waves) 

Hydrodynamic shallow water waves are described by Euler's equation: 

dv ^ 1 

— + vS/v = VP. 

ot p 

The fluid pressure is 

P^Po + P9{h-z), 

where Pq is constant, such that the x- and y-components of the Euler equations are 

dh 
^ dx^ 
dh 
dy 

An additional equation of motion comes from the incompressibility condition: 



dvx 


dvx 
dx 


dvx 
dy 


dt 


dvy 


dvy 
dx 


dvy 
dy 


dt 



(Al) 
(A2) 

(A3) 



MHD Shallow Water Waves 

But we also have 



Dh dh dh dh , 



It follows that 

dh dh dh / dv^ dvy 
dt ^ dx ^ dy \ dx dy 



or, more compactly, 

|^+V.(H-0. (A7) 



dt 


dx 


dy 


dh' 


= 0, 


dt 


dx 


dy 


dh' 


= 0, 


dh' 


C)h' 

dx 


V 4 
^ dy 


(dv'. 




'dt 


\ dx 


dy J 



One can easily check that the derived equations in (|A3p and (|A6p satisfy "state of calm" conditions, 

v, = Vy^Q,h = Ho, (A8) 

where Hq is a constant. 

Considering small perturbations to the velocity and height of the water waves yield the following linearized equations: 



(A9) 

H ^Q. 

° dy \ (Jx cly / 
Seeking wave solutions to equation (|A9p yields: 

v,„A(3 + hagk, = 0, 

VygAo + hogky = 0, (AlO) 
v,„k,H + VygkyH + hoAo = 0. 

The linear, simultaneous equations in (jAlOp can be arranged in the form (i.e., "Cramer's Rule"; e.g., Arfken & 
Weber 1995): 

A\ Vy„ =0, 

where 

/ Ao gk, \ 
A = [ Ao gky ] . 
V k^H kyH Ao / 

The solutions for v^a , Vy^ and h^ are non-trivial only if 

dcti = 0. (All) 

It follows that 

Ao [Al - gH {kl + kl)] = 0, (A12) 
which yields the following dispersion relations: 

' (A13) 



\k,V, + kyVy ± ^gH {kl + kl). 

The first of these dispersion relations is trivial; it simply describes an advection wave. Setting Vx — Vy = and 
analyzing the w > solutions, we obtain the x- and y-components of the phase velocity: 



Cp. = ^^/kx = Co Jl + fc^/fc2, 

\ (A14) 

Cpy = uj/ky = coy^l + 



such that 



Cp = + cl - co^2 + kl/kl + kl/kl. (A15) 
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Similarly, the x- and y-components of the group velocity are: 



du! 


kxCo 




dkx 




kl 


du! 


kyCQ 




dky 




— 1 

1,2 

Ky 



(A16) 



and 



^9 = V^L =^0- (^17) 

Note that Cg < Cp and the amplitudes w^o, Vyo ^-nd ho are always in phase. 

A. 2. With Constant Rotation: f -Plane Treatment (Poincare Waves) 
In the case of non-zero rotation, the Euler equation in the rotating frame has a Coriolis force term: 

^ +u.Vt7= -iyP- 2 (f1 X i/) . (A18) 



dt p 
Assuming that 51 (0, 0, f2), the equations of motion become 



dvx 


dvx 


dvx 


dh 


dt 


dx 


dy 


-'Tx 


dvy 


dVy 


dvy 


dh 


dt 


^^^^ ' 

dx 


dy 


dy 



(A19) 

with equation (|A6|) remaining unchanged. Such a treatment with the Coriolis parameter, / — 2f2, held constant 
is called the "/-plane approximation" (see G82 and KC04). It is possible to Taylor-expand / and keep the first 
order term, which will have a spatial dependence. Such a treatment is known as the "/3-plane approximation" and is 
examined in Appendix IA.3I 
We again consider small perturbations to the physical quantities and linearize the equations: 



(A20) 



dt 


Wx 


dx 


\-Vy 


dy 


dh' „^ , 


dv'y 

dt 


Wx 


^. 

dx 




dv^ 
dy 


dh' 

f,^ + 2fl< 



where the linearized equation of motion for h' remains unchanged (see equation |A9] ) . 
The matrix A now becomes 

CiAo —2fl igkx \ 
20. iAQ igky I , 
kxH kyH Ao J 

and the non-trivial dispersion relation is 

w = kxVx + kyVy ± ^gH{kl + kl)+An-^. (A21) 

Such a dispersion relation agrees with that described by Holton (2004; pg. 208). Waves with this dispersion relation 
are often called "Poincare waves" (pg. 196 of G82), while Poincare waves with a finite horizontal boundary are termed 
"Kelvin waves" (pg. 615 ofKC04)6. 
Setting Vx = Vy = (i.e., solid body rotation), the phase velocity is 

Cp = ^gH (2 + kl/kl + kl/kl) + 4fl2 {i/kl + l/fc2), (A22) 

while the group velocity is 



Co 



fc2 



k^ + AVL'^/gH 



1/2 



(A23) 



Poincare waves are often referred to as "gravity waves" or "g-modes" with rotational effects implied. 
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A. 3. With Constant Rotation: fi-Plane Treatment (Poincare and Rossby Waves) 
Consider the Coriolis parameter to be 

/ = /o + /3y, (A24) 
where /o = 2f7 sin 8 and Q denotes the latitude on a sphere. The parameter f3 is given by 

/3=— cose, (A25) 
K 

with R being the radius of the sphere. Such a treatment is known as the "/3-plane approximation" (e.g., KC04), and 
is the first step towards considering curvature effects on a sphere. 
The set of linearized equations in the /3-plane approximation is: 

5< ^ dh' 

dh' , , „ 

^+5^ + K = 0, (A26) 
dh' fdv'^ dv'y 



dt ' \ dx ' dy ' ^ 
Differentiating the second equation in (|A26|) twice with respect to t yields: 



dt^ dt dx \dx dy J dydt \ dx dy 

Differentiating the first and second equations in (jA26|) with respect to y and x respectively, and subtracting them from 
one another gives: 

d fdv' dv'\ Jdv' dv' 



Adding equations l|A27p and (jA28|) . with the latter operated on by gH-^, eliminates v'^: 

d^v' ,dv' d f d^v' d^v'\ dv' ^ ^ 

-w^^^- -^^^ + # j - = 

Such an approach is described in KC04 (pg. 610); a generalized version of it was applied to MHD systems by 
Zaqarashvili et al. (2007). From equation (IA29|1 . one can seek wave solutions as usual and obtain the dispersion 
relation. 

There is a more convenient approach one can adopt to obtain the dispersion relation. Instead of juggling differential 
operators, one differentiates the first equation in (jA26p with respect to y; the purpose of this step is to extract the 
first order terms in /. Together with the other two equations in (jA26|) . wave solutions are then sought and the matrix 

A is once again constructed: 

/ ujky - (/3 + iky fa) -gk^ky \ 
A = I fo -iLu igky 



kxH kyH ~uj 



Setting detA — yields the dispersion relation: 



UJ 



3 , , I „TT1„2 



{gHe + /2) - gHk^p = 0. (A30) 

Notice that equation (|A30p is asymmetric with respect to k^ and fcj,, implying that wave motion is not isotropic in 
the horizontal direction. Very slow (w ^ /), large-wavelength {kH <^ 1) waves are called "Rossby waves". They are 
sometimes called "planetary waves" , and have an angular frequency of (see also pg. 634 of KC04) 

2- (A31) 



fc^ + (/o/co) 

For very large wavelengths, the eastward phase speed is 



CR«-/3(co//o)'. (A32) 



The negative sign shows that phase propagation is always westward. 

B. CLASSICAL SHALLOW WATER WAVES: SPHERICAL COORDINATES 
We revisit the classic work of LH68 as a prelude to [jS] 
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B.l. Equations 

Consider waves with a velocity v — {vr, ve, v^). The hnearized equations of motion have the form: 

dv'n „^ , „ a dh' 
^ - 217u^cos6l + f;— = 0, 



dt 

+ 2VLv'g COS 

dh' 
'dt 



H 



RsinO 



R 89 

g dh' 
Rs'mO d(f) 



0, 



(Bl) 



= 0. 



The amphtude equations become: 



mveo + HV4,g + DriQ = 0, 
IJ.ve„ + wv^g - sryo = 0, 

zue (l - n^) rjo - Dvg^ - sv^^ 



(B2) 



0. 



The second equation in the preceding set can be used to ehminate and obtain LgVgg — 0, where the operator is 



2e/itx7 ( wD — sfj, 



hero 



(1 - ^2) 



(B3) 



in agreement with equation (7.8) of LH68 (see Appendix [C|) . 



B.2. Near-Equator Solutions (e^ 1) 

When e ^ 1, the last term in equation (|B3P is small compared to the other terms. The first and second terms cannot 
be dropped because the boundary conditions at /i = ±1 have to be satisfied (LH68)^. In this approximation, one may 

recognize LsVeo = as the "spheroidal wave equation" with the solutions Ssn (v^jM) (Abramowitz & Stegun 1970), 
where the "separation constant" is 



As„ (g) = ero^ 



and q = ^/e. Note that for I = n ~ s, 



^sn (q) 



21 + 1 



(B4) 



(B5) 



for g > 4 (Abramowitz & Stegun 1970). When e = 0, Lg becomes the operator for the associated Legendre equation; 

the corresponding solution is Vn\tj) and the separation constant is As,i — n{n + 1) (Arfken & Weber 1995). In 
this case, only the eastward-propagating Poincare mode exists. When e 7^ 0, a westward-propagating Poincare mode 
appears in addition to the Rossby mode. The dispersion relation is extracted from equation (|B4[) : 



njAsri ( 



.1^1 =0. 



(B6) 



The Poincare and Rossby waves are described by the so-called "Type 1" and "Type 2" solutions of LH68, respectively. 
The azimuthal velocity perturbation in both cases has the functional form, 

r 



where 7ij(X) is the Hermite polynomial and 



= exp 

We first focus on the Type 1 solutions, where 



2 



UJ ^ i: 



exp i [scf) — Lot). 



Making use of the second and third equations in the set (jB2p leads to 



ru'^e (1 — fj?) — 



(B7) 



(B8) 



(B9) 



(BIO) 



Specifically, the spheroidal wave functions need to be finite at the poles (Abramowitz & Stegun 1970). 
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The preceding approximation is justified because siivg^ <ti wDvg^,, and w'^e cx e^/^ ^ s^. As we will see later, these 
approximations are inviable for the Type 2 solutions. It follows that 

Notice that we have the "±" instead of the sign in equation (8.32) of LH68, because we are defining as the 
co-latitude. 

To obtain v'^, we use the second equation in (jB2p . 

-.0 - (B12) 
w 

The preceding approximation is valid because /iweg cx: e^^l^ ^ srj cx e~^/'^. It follows that 

'^^ ^mi-i + \ui+^ (B13) 



Obtaining the Type 2 solutions, where 
involves realizing that the approximations made to obtain the Type 1 solutions do not hold, because 

-1/4 rS~ -1/4 

sfiVgg CX e ' ^ ruUvgg cx e ' , 

(B15) 



2 2 



U7 W 



and therefore the previously dropped terms need to be retained. Making use of the recursion relations for Hermite 
polynomials, we arrive at: 



. 2^ + 1 1 



(B16) 



Again, notice that the plus and minus signs are switched compared to LH68. Note that we require I > and ^ > 1 for 
Type 1 and 2 solutions, respectively. 

B.3. Near-Pole Solutions f— e ^ 1) 
LH68 showed that ru — > ±1 as —e increases. He reasoned that for the penultimate term in equation (|B3p to not 
dominate the other terms, we must have — e (tu^ — /i^) ~ 1. When — e ^ 1, (1 — ^P) <C 1 and /i^ ^ 1. Therefore, the 

solutions to LgVOg = are concentrated near the poles for large values of — e. 
To proceed, one assumes va to have the form, 

vj = ±1 + (B17) 



Using the substitution, 

Y = V~e (1 - fi^) , (B18) 



the equation LsVeg = becomes 

1 Q s^±2s 



_dF2 4 ' 2Y 4^2 
readily recognized as Whittaker's equation (Abramowitz & Stegun 1970). The solutions are 



ieo = 0, (B19) 



Vo exp (^-^^ y(™+i)/2 (y) ^ (B20) 



where c\™\y) is the associated Laguerre polynomial (Arfken & Weber 1995). For m w ±1, we have m = |s ± 1|; also, 
Q = 2m + 2^ + 1. The n7 « — 1 and +1 cases correspond to the "Type 4" and "Type 5" so lutio ns of LH68, respectively. 
The quantities h' and v'^ can be calculated using the procedure described in Appendix IB.2I to compute the Type 2 
solutions. 
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C. GENERALIZED IDENTITY OF LONGUET-HIGGINS (1968) 

We provide a proof for the identity used to go from equation (7.5) to (7.8) of LH68 and generalize it for MHD 
systems. Using the test functions T = .?^(m) ^'^^ S — GifJ-), one show can that 



(x-D + Ms) [g (xD - ms)] ^ = a (x^ + Ms) (xD -^ls^T + x (dg) [x (dt) - fisT 
(xD + ^is) (xD - ^ - X^D (dt) - xD (fisT) + ^isx (dt) - {fisf T. 



Now, we consider 
such that 



G = [n7xe (1 - 
2wx^^i 



xD-^is]T + uJe (x'-m')-^ = 0. 



tuxe (1 — /i^) — 

The operator used in Longuet-Higgins (1965) and LH68, originally denoted by V^, is defined by us as 



1 -^^ 



(CI) 
(C2) 

(C3) 
(C4) 

(C5) 



To use these identities for non-magnetized systems, simply set x = ^• 

D. /3-PLANE TREATMENT FOR MHD EQUATORIAL EIGENFUNCTIONS 

The equatorial eigenfunctions can be obtained in the /3-plane approximation, an analysis first performed by Matsuno 
(1966), who recognized it as an eigenvalue problem. The /3-plane analysis is algebraically more tractable and provides 
a consistency check on the results in ij3.2l We first normalize the equations using the length and time scales: 

L= ,fM., T=^^. (Dl) 



13' ' V^' 

We also let 70 = gho/c^^^p-^^^, where ho = ho/L. Subsequently, quantities in this section marked with a "tilde" are 
normalized by L and/or T, unless otherwise specified. 

Seeking expilk^x — cot) wave solutions, the amplitude equations become: 

dy 

iXVxo - yvyo + ikxlQ = O7 (D2) 



0, 



dy 



1(1)70 = 0, 



where 



X = ^ 

= 



(D3) 



Manipulating the second and third equations in (|D2p yields: 

1 f~ dvy^ 



70 



I {ibx - ~ki) V dy 



dy 
dv 



■ i^yVyi 



(D4) 



Substituting the preceding expressions into the first equation of (jD2p gives: 

d^ 



dY"^ 



'1/2 



0. 



(D5) 



As expected, the governing equation for Vy^ is the parabolic cylinder equation, where a = w/x and Y = a^^'^y. It has 
the solution: 



= exp 



Vo Hi 



expi(fca;X — Lot). 



(D6) 
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Here, TLi = TLiiY) and Vq is an arbitrary normalization constant. 
The dispersion relation is: 



6Bi{Bi + kt +kt-{2l + iy 



-2Bl (2Bt 



kUl-6Bt +Bii2l + iy 



(D7) 



+ 2k^ {2Bl + kl^ Co^ - 2hBl {bI + kl 



U! 



Bt { Bt + 2Bfkl + kt) =0 



For the magneto-Rossby and west magneto-Poincare modes, the other eigenfunction solutions are: 



l' 



i [ujx - 



ini-iB+ - 



(D8) 



where 7' = gh! j t?^"^ 0^1"^ . For the east magneto-Poincare mode, we have: 

^0 



7 



Analogous to the case of spherical geometry, we have: 



IHi-iA- 



-H1+1B+ 



-'H1+1A+ 



(D9) 



A± = k, 
B± = kx 



1/4 



1/4 



(DIO) 



±1x1 



r.|V4|~|3/4 



The (normalized) magnetic field perturbations are obtained fi'om hx,y — iBzv',^ yl^- 
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